function [LL, varargout] = loglikelihood_memory_efficient(modelId, data, params, takeNegative)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% A version of loglikelihood that is less demanding on memory that needs to be allocated at any given time
	% (when combined with loglikelihood_by_blocks).
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% modelId:			       integer (1 => BinomialLogistic, 2 => Poisson)
	% data:			           object:
	% 	.dims:				     	object:
	%		.dims1:						cell(NumParts,1)
	%			{ii}:						integer: gives dim1_ii that corresponds to Xparts{ii}
	%		.Xpart_2_NumX:				1 x NumParts: gives integers (dim2 of Xparts{ii}.X)
	%		.Xpart_2_NumX_FEs:			1 x NumParts: gives integers (dim2 of Xparts{ii}.X_FEs)
	%		.Xpart_2_Num_FE_vals:		cell(1,NumParts)
	%			{ii}:						1 x NumX_FEs_i  --> gives integer: number of possible values for Xparts{ii}.X_FEs(:,ff)
	%		.NumFEvals2Keep:			cell(1,NumParts) --> gives integer
	%		.NumParts
	%		.NumObs
	%		.NumParams
	%    .LL_offset:               1 x 1
	%	 .Xparts:				   cell(NumParts,1)
	%	 	{ii}:				       object
	%			.X:					       dim1_ii x NumX
	%			.X_FEs:					   dim1_ii x NumX_FEs
	%			.NumX_FE_vals:			   1 x NumX_FEs: gives integer
	%    .N:                       NumObs x 1    (if modelId==1)
	%    .logLambdaOffset_dim1:         [T x K1]      (if modelId==2 or modelId==3)
	%    .logLambdaOffset_dim2:     [T x 1 x K2]  (if modelId==2 or modelId==3)
	%    .Y:                       NumObs x 1
	%    .spellDurations:          T x 1         (if modelId==3)
	% params:                  NumX x 1
	% takeNegative:            boolean
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% LL:            		scalar
	% grad:					NumParams x 1
	% hessian:				NumParams x NumParams
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	dims = data.dims;
	
	if (modelId == 2 || modelId == 3) && dims.dimType == 2
		K = size(data.logLambdaOffset, 2);
		logLambdaOffset = data.logLambdaOffset; % T x K
		logLambdaOffset = reshape(logLambdaOffset, [prod(size(logLambdaOffset)) 1]); % (T*K) x 1
	end
	
	if (modelId == 2 || modelId == 3) && dims.dimType == 3
		logLambdaOffset = data.logLambdaOffset_dim1 + data.logLambdaOffset_dim2; % T x K1 x K2
		logLambdaOffset = reshape(logLambdaOffset, [prod(size(logLambdaOffset)) 1]); % (T*K1*K2) x 1
	end
	
	if modelId == 3 && dims.dimType == 2
		spellDurations = repmat(data.spellDurations, [1, K]); % T x K
		spellDurations = spellDurations(:); % (T*K) x 1
	end
	
	if modelId == 3 && dims.dimType == 3
		K1K2 = dims.dims1{1};
		spellDurations = repmat(data.spellDurations, [1, K1K2]); % T x (K1*K2)
		spellDurations = spellDurations(:); % (T*K1*K2) x 1
	end
	
	% Split params into parts
	params_parts = split_combine_parts(dims, 1, params);
	
	% Compute V
	V = compute_Xbeta2(dims, data.Xparts, params_parts);	% NumObs x 1
	
	%%% Compute LL, dLL_dV and d2LL_dV2 (model specific)
	if modelId == 1  % Logistic-Binomial
		error('This function has not been written for modelId=1.');
	end
	if modelId == 2  % Poisson
		logLambda  = logLambdaOffset + V;								% NumObs x 1
		clear V logLambdaOffset;
		lambda     = exp(logLambda);										% NumObs x 1
		LL         = data.LL_offset  + data.Y'*logLambda -sum(lambda,1);	% 1 x 1
		clear logLambda;
	end
	if modelId == 3  % Exponential
		logLambda  = logLambdaOffset + V;						% NumObs x 1
		clear V logLambdaOffset;
		D_lambda = spellDurations.*exp(logLambda);				% NumObs x 1
		clear spellDurations;
		LL         = data.LL_offset  + data.Y'*logLambda - sum(D_lambda);	% 1 x 1
		clear logLambda;
	end
	
	%%% Output LL and derivatives
	if takeNegative; LL = -LL; end;
	
	if nargout >= 2
		%% Compute dLL_dV
		if modelId == 1  % Logistic-Binomial
			error('This function has not been written for modelId=1.');
		end
		if modelId == 2  % Poisson
			dLL_dV = data.Y - lambda; 			% NumObs x 1
		end
		if modelId == 3  % Exponential
			dLL_dV = data.Y - D_lambda; 			% NumObs x 1
		end

		%% Compute grad
		grad = compute_XprimeY2(dims, data.Xparts, dLL_dV)';		% NumParams x 1
		clear dLL_dV;
		if takeNegative; grad = -grad; end;
		varargout{1} = grad;
	end
	
	if nargout >= 3
		%% Compute d2LL_dV2
		if modelId == 1  % Logistic-Binomial
			error('This function has not been written for modelId=1.');
		end
		if modelId == 2  % Poisson
			d2LL_dV2 = -lambda; % NumObs x 1
		end
		if modelId == 3  % Exponential
			d2LL_dV2 = -D_lambda; % NumObs x 1
		end
		
		hessian = compute_XprimeY_X2(dims, data.Xparts, d2LL_dV2);	% NumParams x NumParams
		clear d2LL_dV2;
		if takeNegative; hessian = -hessian; end;
		varargout{2} = hessian;
	end
end
